SUBROUTINE HESS (NPARS,PARVEC,H,LDH)
USE SHARE
IMPLICIT NONE
INTEGER :: i,j,ja,jb,ip,ia,LDH,NPARS,jp
DOUBLE PRECISION, DIMENSION(900,10) :: num
DOUBLE PRECISION, DIMENSION(900) :: den
DOUBLE PRECISION :: Den_hes(900),mat(10,10,900),wt(6000)
DOUBLE PRECISION :: PARVEC(10)
DOUBLE PRECISION, DIMENSION (LDH,LDH) :: h
! PRINT 7685, LDH,NPARS; 7685 format (' in hess ldh ',i4,' npars ',i4,' NMOS ',i4)
DO i=1,NMOS   !  compute the summation parts of denominator
   den(i)=zero
   DO jp=1,NPARS
      num(i,jp)=zero
   END DO
   DO j=1, ktr(i)
         DO jp=1,NPARS
            num(i,jp)=num(i,jp)+chars(j,i,jp)*rit(j,i)
         END DO
         sc=zero
         DO jp=1,NPARS
            sc=sc+(parvec(jp)*chars(j,i,jp))
         END DO
         wt(j)= weig(j,i)+ sc/ktr(i)
         den(i)=den(i)+wt(j)*rit(j,i)
   END DO
   DO jp=1,NPARS
      num(i,jp)=num(i,jp)/ktr(i)
   END DO
END DO
DO i=1,NMOS  !  denominator for gradient and hessian in each month
   Den_hes(i) = (1.0d0+den(i))**(-gamma-one)
END DO
DO i =1,NMOS  ! gradient and hessian in each month
   DO ip=1,NPARS
      DO ja=1,NPARS
         mat(ip,ja,i)=gamma*Den_hes(i)*num(i,ip)*num(i,ja)/NMOS
      END DO
   END DO
END DO
DO ip=1,NPARS
   DO jp=1,NPARS
      h(ip,jp)=zero
   END DO
END DO
DO i=1,NMOS
   DO ia=1,NPARS
      DO ja=1,NPARS
         h(ia,ja)=h(ia,ja)+mat(ia,ja,i)
      END DO
   END DO
END DO
RETURN; END
